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ABSTRACT Bovine nnastitis is an inflannnnation-driven disease of the bovine nnannnnary gland that costs the 
global dairy industry several billion dollars per year. Because disease susceptibility is a multifactorial complex 
phenotype, an integrative biology approach is required to dissect the molecular networks involved. Here, we 
report such an approach using next-generation sequencing combined with advanced network and pathway 
biology methods to simultaneously profile mRNA and miRNA expression at multiple time points (0, 1 2, 24, 36 and 
48 hr) in milk and blood FACS-isolated CD14+ monocytes from animals infected in v/Vo with Streptococcus uberis. 
More than 3700 differentially expressed (DE) genes were identified in milk-isolated monocytes (MIMs), a key 
immune cell recruited to the site of infection during mastitis. Upregulated genes were significantly enriched for 
inflammatory pathways, whereas downregulated genes were enriched for nonglycolytic metabolic pathways. 
Monocyte transcriptional changes in the blood, however, were more subtle but highlighted the impact of this 
infection systemically. Genes upregulated in blood-isolated monocytes (BIMs) showed a significant association 
with interferon and chemokine signaling. Furthermore, 26 miRNAs were DE in MIMs and three were DE in BIMs. 
Pathway analysis revealed that predicted targets of downregulated miRNAs were highly enriched for roles in 
innate immunity (FDR < 3.4E-8), particularly TLR signaling, whereas upregulated miRNAs preferentially targeted 
genes involved in metabolism. We conclude that during S. uberis infection miRNAs are key amplifiers of mono- 
cyte inflammatory response networks and repressors of several metabolic pathways. 
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Bovine mastitis is an inflammation-driven disease of the bovine mam- 
mary gland that is associated with significant costs to the global dairy 
industry. In Europe this cost is estimated to be approximately €2 
billion euro per year (Wells et al. 1998), with similar figures available 
in the United States (Jones and Bailey 2009). Causative agents of mastitis 
infection include, but are not limited to, coliforms (Escherichia coli), 
Streptococci (Streptococcus uberis) and Staphylococci (Staphylococcus 
aureus). S. uberis is now ranked among the most prevalent mastitis- 
causing pathogens throughout Europe and North America (Reinoso 
et al. 2011; Ward et al. 2009). 

Mastitis develops as bacteria entering the udder via the teat canal 
stimulate a pathological form of inflammation. Bacteria encounter 
epithelial cells lining the mammary gland stimulating a local inflammatory 



;-2iG3*Genes | Genomes | Genetics 



Volume 4 I June 2014 I 957 



response that facilitates their transport across the epithehal barrier, 
where they are detected by resident immune cells, such as monocytes. 
Both cell types constitutively express surface pathogen recognition 
molecules such as Toll-like receptors (TLRs), enabling them to 
function in a sentinel capacity. Invasive S. uberis triggers TLR2 and 
TLR4 mobilizing local and systemic inflammatory mediators 
(Bannerman et al. 2004; Moyes et al. 2009). Typically, chemokines, 
interleukins (ILs), and tumor necrosis factor (TNF)-a initiate local 
physiological changes in vascular permeability, cell differentiation, 
and apoptosis. Concurrently, systemic innate immune changes 
provoke acute phase protein (APP) production, which is distributed 
systemically to suppress the spread of bacteria locally (Mitterhuemer 
et al. 2010). During this phase, immune cells are recruited to the point 
of infection (Rinaldi etal 2010). Monocytes are released from the bone 
marrow into the circulatory system and eventually reach the 
mammary gland via chemokine ligand-mediated cell migration. 
There they differentiate into macrophage and dendritic cell popula- 
tions (Dong et al. 2013; Shi and Pamer 2011). Neutrophils comprise 
the majority of immune cells in an infected gland during an infection. 
Neutrophils are tasked with directly clearing invasive bacteria via 
phagocytosis or neutrophil extracellular traps (NETS) and subse- 
quently aid in resolution of inflammation (Lippolis et al. 2006; 
Reinhardt et al. 2013). Once recruited to the site of infection, mono- 
cytes and neutrophils orchestrate antimicrobial activity to control 
bacterial spread and resolve the infection (Dong et al. 2013; Serbina 
et al. 2008). Immune cells and other somatic cells can be detected in the 
milk of infected animals, and the counts of the number of such somatic 
cells per ml, called the somatic cell count, is an indicator of mastitis 
(Jones and Bailey 2009). 

The local immune response in mammary tissues has been 
examined by several approaches both in vivo and in vitro. Candidate 
gene-based approaches and microarray technology have determined 
that more than 2000 genes spanning immunity, metabolism, and 
tissue remodeling are active during mastitis (Mitterhuemer et al. 
2010; Moyes et al. 2009; Swanson et al. 2009). However, modest data 
are available examining transcriptional activity in either milk or blood 
monocytes from infected animals (Prgomet et al. 2005), and little is 
known regarding the role microRNAs (miRNAs) play in regulating 
these responses. 

The miRNAs are small, noncoding RNAs that play a key role in 
the regulation of innate and adaptive immunity as posttranscriptional 
regulators of gene expression (O'Connell et al. 2010). They have been 
shown to regulate immune function in several cell types. Neutrophil 
senescence, for example, is regulated by a discrete miRNA repertoire 
(Ward et al. 2011). Naive mouse B cells are indirectly regulated by 
miR-155 via histone deacetylase 4 repression, whereas naive CD4+ T- 
cell differentiation and function are regulated by global changes in 
miRNAs (Bronevetsky et al. 2013; Sandhu et al. 2012). In a recent 
study, we concluded that miRNAs likely play a key role in regulating 
the innate immune response in mammary epithelial cells to a bovine 
mastitis pathogen in vitro (Lawless et al. 2013). 

Although miRNA expression is abundant in numerous bovine 
tissues, genome-wide studies elucidating the regulatory roles of 
miRNAs in bovine immunity are limited (Coutinho et al. 2007; Jin 
et al. 2009; Xu et al. 2009). Furthermore, no bovine studies to date 
have applied next-generation sequencing (NGS) to examine global 
miRNA expression in immunity and infection in vivo. In this study, 
we report a NGS approach to profile the expression of bovine 
miRNAs and mRNAs at multiple time points in milk and blood 
CD 14+ monocyte cells isolated from Holstein Friesians infected 
in vivo with S. uberis, a causative agent of bovine mastitis. 



MATERIALS AND METHODS 
Animals 

Ten female Holstein Friesians in the middle of their first lactation 
period, aged between 26 and 30 months and between 3 and 5 months 
postpartum, were selected for this study. The trial was conducted at 
the USDA National Animal Disease Center (NADC) in Ames, Iowa. 
AU animals had a medical history that was free from mastitis. The 
National Animal Disease Center's Animal Care and Use Committee 
approved all procedures used in this study. 

Infection protocol 

Five animals were infected via the teat canal of the right front quarter 
with approximately 500 colony-forming units (CFU) of a mastitis- 
causing pathogen, S. uberis 0140, in 10 ml saline. Five control animals 
were inoculated with saline only. Milk and blood samples were 
obtained from each animal at 0, 12, 24, 36, and 48 hr after infection 
(or mock infection) as described below. At each time point, rectal 
temperature, total volume of milk, somatic cell count, bacterial counts, 
ambient temperature, humidity, and additional observations were 
recorded for each animal. Bacterial counts were determined from 
5 -ml milk samples collected aseptically from the infected quarter. Milk 
was serially diluted in sterile phosphate-buffered saline and spread on 
blood agar plates and then incubated for 24 hr at 37°. After incuba- 
tion, plates were examined for bacterial growth and CPUs per ml were 
determined. 

Cell extraction from milk 

Milk was collected using a sterilized quarter milker from infected and 
control animals at each time point. The total volume of milk was 
noted; 5 ml milk was isolated for milk bacteriology. The remaining 
milk was then diluted into Hanks' Balanced Salt Solution (HBSS) 
without Phenol red, Mg, and Ca... plus 10 mM EDTA. The mixture 
was inverted several times, transferred to a 1 -liter centrifuge bottle, 
and centrifuged in a fixed angle rotor at 10,000 x g for 30 min. After 
centrifugation, the supernatant was poured off and the pellets were 
resuspended in 150 ml HBSS plus 5 mM EDTA. The resuspended 
pellets were then transferred to fresh centrifuge tubes and spun at 
2500 X g for 30 min. After the second spin, the supernatant was 
poured off and the pellet from each sample was resuspended in 20 
ml RPMI 1640 plus 1 mM sodium pyruvate plus 2 mM L-glutamine 
plus 50 ug/ml gentamycin plus 10% FBS (cRPMI) (Sigma- Aldrich, 
Steinheim, Germany). The 20-ml cell suspension was divided into 
2 X 10 ml aliquots in 15 ml conical centrifuge tubes. The cells were 
pelleted by centrifuge at 650 x g. After this spin, the pellets were 
pooled into 4 ml cRPMI and were labeled for cell sorting. The cells 
were counted by trypsin blue exclusion (Careforde, Chicago, IL) to 
determine the cell count. 

Cell extraction from blood 

At each time point, animals were led into a crush, and 2x60 ml 
syringes of blood were extracted by venipuncture and immediately 
placed on ice. The total volume of blood and total number of cells 
in blood were determined using a hemocytometer. Blood was spun for 
20 min at 1200 X g. The buffy coat was observed between serum and 
red blood cell phase and was removed. Contaminating red blood cells 
were lysed by adding 1 volume lysis solution (10.6 mM Na2HP04 and 
2.7 mM NaH2P04) and inverting the tubes several times, immediately 
followed by adding half volume restore solution (10.6 mM Na2HP04, 
2.7 mM NaH2P04, and 460 mM NaCl) and inverting the tubes. The 
buffy coat was then spun for 10 min at 650 x g, and the red 
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supernatant was poured off. Cells were resuspended in red blood cell 
lysis solution and then restore solution again as described above. Cell 
were spun for another 5 min at 650 X g and resuspended in 5-10 ml 
media. 

Isolation of CD14+ monocytes by flow cytometry 

Milk-derived and blood-derived CD 14+ monocytes were isolated by 
fluorescence-activated cell sorting (FACS). Briefly, cells were labeled 
with monoclonal anti-bovine CD14 (Clone CAM36A; VMRD, Pullman, 
WA) and a PE-conjugated anti-mouse IgGl antibody (Southern 
Biotechnology, Birmingham, AL). Labeled cells were separated based 
on fluorescence intensity using the BD FACS Aria Cell Sorting System 
(BD Biosciences, San Jose, CA). Cells with more than 95% purity were 
isolated from the milk and peripheral blood of each animal. 

mRNA extraction 

The mirVana RNA Isolation Kit (Ambion, Austin, TX) was used to 
extract total RNA from FACS-isolated cell populations. Procedures 
were performed according to the manufacturer's protocol (Supporting 
Information, File SI). RNA was quantified and integrity was con- 
firmed using an Agilent RNA Kit on a 2100 Bioanalyzer platform 
(Agilent Technologies, Loveland, CO) (File SI). 

miRNA extraction 

MicroRNA was extracted using mirPremier microRNA Isolation Kits 
( Sigma- Aldrich, Steinheim, Germany). Procedures were performed 
according to the manufacturer's protocol (File SI). Small RNA was 
quantified using an Agilent small RNA Kit on a 2100 Bioanalyzer 
platform (File SI). 

mRNA library generation 

One hundred indexed mRNA libraries (50 blood monocyte libraries 
and 50 milk monocyte libraries) were prepared for cluster generation 
using TruSeq v2 RNA sample preparation kits (Illumina, San Diego, 
CA). Procedures were performed according to the manufacturer's 
protocol (File SI). The finished libraries were validated on an Agilent 
Bioanalyzer 2100 using an Agilent DNA-1000 chip (Agilent, Loveland, 
CO), at which point they were loaded for cluster generation. The 
samples were sequenced on an Illumina HiSequation 2000 at the Iowa 
State Sequencing Center (50 -bp single- end). Infected and control sam- 
ples (n = 100) were randomized across four flow cells (/.e., three or 
four samples multiplexed per lane) to avoid confounding flow cell/ 
lane effects (Auer and Doerge 2010). The barcode compatibility chart 
provided with the TruSeq RNA sample preparation kit was adhered to 
when pooling libraries. Fastq files were produced using the CASAVA 
1.8 pipeline. 

mRNAseq analysis 

One hundred fastq files were generated containing the sequencing 
data for each of the 100 mRNAseq libraries. The quality and number 
of the reads for each sample were then assessed using FASTQC 
vO. 10.0 (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/). Each 
sample was then put through a number of quality- control filters. First, 
reads were filtered using the fastq Illumina filter vO.l (http://cancan. 
cshl.edu/labmembers/gordon/fastq_illumina_filter/), which removes 
reads from the fastq files that were flagged as not passing the lUumina 
CASAVA pipeline filters. Cutadapt vl.2 (http://code.google.eom/p/ 
cutadapt/) was used to trim the adaptors from reads when necessary. 
The remaining reads were then further filtered using the fastq qual- 
ity filter package (http://hannonlab.cshl.edu/fastx_toolkit/) vO.0.13.2. 



When at least 70% of the bases had a Phred score <20, reads were 
removed. Reads passing all the aforementioned filters were also 
trimmed at their ends to remove low-quality bases (Phred score 
<20). Reads that were <20 nt after trimming were discarded. Reads 
that passed all quality- control steps were then aligned to the bovine 
genome (UMD3.1 assembly) (Zimin et al 2009) using TopHat v2.0.8 
(Trapnell et al 2009), allowing one mismatch. Reads that did not 
uniquely align to the genome were discarded. HTSeq- count version 
0 .5 .3p3 (http :/ / wwwhuber. embl.de/ users/ anders/HTSeq/ doc/ overview, 
html) using the union model was used to assign uniquely aligned 
reads to Ensembl (v69) annotated bovine genes. 

Differential gene expression analysis 

Data were normalized across libraries using the trimmed mean of M- 
values (TMM) normalization method (Bullard et al 2010). The R 
(version 2.15.2) Bioconductor package EdgeR (v2.4.6) (Robinson 
et al 2010), which uses a negative binomial distribution model to 
account for both biological and technical variability, was applied to 
identify statistically significant differentially expressed (DE) genes. 
Any samples that had less than five million uniquely aligning reads 
were removed from further analysis. Only genes that had at least one 
count per million in at least three samples were analyzed for evidence 
of differential gene expression. The analysis was undertaken using 
moderated tagwise dispersions. DE genes were defined as having a fold 
change in gene expression >1.5 and a Benjamini and Hochberg cor- 
rected FDR of <0.05 (Benjamini and Hochberg 1995). 

Hierarchical clustering 

Hierarchical clustering of milk and blood mRNA normalized read 
counts were performed in the R (version 2.15.2) hclust package. Heat- 
maps were generated using the R heatmap package. 

Gene ontology and pathway analysis 

The R (version 2.15.2) Bioconductor package GOseq (version 1.10.0), 
which corrects for gene length bias (Young et al 2010), was used to 
identify overrepresented pathways using pathway annotation imported 
from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa 
et al 2007) database. KEGG disease pathways were excluded to focus the 
analysis on primary signaling pathways. Pathways were considered 
significantly overrepresented with an FDR <0.05. Pathway analysis 
was undertaken using Ensembl predicted human 1:1 orthologs of the 
bovine DE genes. 

Additionally, we manually generated two pathway annotations 
that were of interest but not annotated in detail in KEGG: the "inflam- 
masome" and "interferon" pathways. Gene IDs for these pathways 
were sourced from SA Biosciences (Qiagen) RT^ Profiler PGR Array 
Human Interferon and Receptors (PAHS-064A) and RT^ Profiler 
PGR Array Human Inflammasome (PAHS-097A) annotations. The 
interferon pathway consisted of 84 genes (Table SI) with expression 
controlled by or involved in cell signaling mediated by interferon 
ligands and receptors, whereas the inflammasome pathway consisted 
of 95 key genes (Table SI) involved in the function of inflammasomes, 
protein complexes involved in innate immunity, and general NOD- 
like receptor (NLR) signaling. 

Network analysis methods 

To generate molecular interaction networks, the human 1:1 orthologs 
of bovine genes that were DE during at least one of the four time 
points in milk-isolated monocytes (MIMs) from S. w^er/5-infected 
animals were uploaded to InnateDB (www.innatedb.com) (Lynn 
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et al. 2008). InnateDB is one of the most comprehensive databases of 
all human and mouse experimentally supported molecular interac- 
tions (> 300,000 interactions in July 2013) but also specifically 
includes annotation on more than 19,000 manually curated human 
and mouse innate immunity relevant interactions, many of which are 
not present in any other database (Lynn et al 2010). Networks were 
visualized using Cytoscape 2.8.2 (Shannon et al. 2003). 

The network was analyzed using the cytoHubba plugin (Lin et al 
2008) for Cytoscape 2.8.2 (Shannon et al 2003) to identify network 
hubs and bottlenecks using default parameters. The jActiveModules 
plugin (Ideker et al 2002) in Cytoscape 2.8.2 (Shannon et al 2003) 
was also used to identify high-scoring DE subnetworks (overlap 
threshold = 0.3; search depth = 3; number of modules = 5; "regional 
scoring" and "adjust score for size" were both enabled). The InnateDB 
pathway analysis tool was used to identify overrepresented pathways 
among module genes. 

miRNA library generation 

One hundred indexed miRNA libraries (50 blood and 50 milk 
monocyte libraries) were also prepared for cluster generation and 
sequencing using the TruSeq Small RNA sample preparation kit. 
These miRNA libraries were prepared from the same samples that the 
mRNAseq libraries were prepared. Procedures were performed 
according to the manufacturer's protocol (File SI). The finished li- 
braries were validated on an Agilent Bioanalyzer 2100 using an Agilent 
DNA high-sensitivity chip (Agilent), at which point they were loaded 
for cluster generation. The samples were sequenced on an Illumina 
HiSequation 2000 (50-bp single-end). Infected and control samples 
(n =100) were randomized across three flow cells (/.e., seven or eight 
samples multiplexed per lane) to avoid confounding flow cell/lane effects 
(Auer and Doerge 2010). Fastq files were produced using the CASAVA 
1.8 pipeline. In a few cases, the sequenced miRNA libraries were found 
to be adaptor-contaminated. These libraries were repurified and rese- 
quenced. The list of resequenced samples can be found in Table S2. 

miRNAseq Analysis 

Preliminary quality- control analysis of the 100 miRNAseq fastq files 
was again performed with FASTQC software vO.10.0 (http://www. 
bioinformatics.babraham.ac.uk/projects/fastqc/). Cutadapt vl.2 
(code.google.eom/p/cutadapt/) was then used to trim 3' adaptor 
sequences. Reads that were shorter than 18 nucleotides after trim- 
ming were discarded. Trimmed reads were then further filtered using 
the fastq quality filter (http://hannonlab.cshl.edu/fastx_toolkit/) 
vO.0.13.2. Reads with at least 70% of the bases with a Phred score 
<20 were removed (Cock et al 2010). Finally, reads passing all these 
filters were also trimmed at their ends to remove low- quality bases 
(Phred score <20). Reads that successfully passed filtering were 
aligned to the bovine genome (UMD3.1) using novoalign v2.08.03 
in miRNA mode (http://www.novocraft.com) allowing one mismatch. 
Nonuniquely aligning reads were discarded. HTSeq version 0.5.3p3 
(http://www-huber.embl.de/ users/ anders/HTSeq/ doc/ overview.html) 
using the union model was used to assign uniquely aligned reads to 
miRBase vl9 miRNA annotation (Kozomara and Griffiths- Jones 
2011). The sequencing data from this pubUcation have been sub- 
mitted to the NCBI GEO database and assigned the identifier 
(GSE51858). 

Differential miRNA expression analysis 

Before assessing differential expression, miRNAseq count data were 
first normalized across libraries using either the trimmed mean of 



M-values (TMM) normalization method (Robinson et al 2010) or the 
upper-quantile normalization (Bullard et al 2010). Differential expres- 
sion analysis of miRNAseq data has been shown to be sensitive to the 
normalization approach implemented (Garmire and Subramaniam 
2012). To address this issue, we identified DE miRNAs in two alter- 
natively normalized datasets: TMM-normalized (Robinson et al 
2010), upper-quantile normalized, and with no normalization. Only 
miRNAs that were identified as DE across all three datasets were 
considered further, i.e., the differential expression of these miRNAs 
was robust to the normalization procedure (Lawless et al 2013). Any 
samples that had less than two million uniquely aligning reads were 
removed from further analysis. The R (version 2.15.2) Bioconductor 
package EdgeR (v2.4.6) (Robinson et al 2010) was applied to identify 
statistically significant DE miRNAs. The analysis was undertaken us- 
ing moderated tagwise dispersions. DE miRNAs were defined as hav- 
ing a Benjamini and Hochberg corrected P value of < 0.05 (Benjamini 
and Hochberg 1995). 

Coexpression and target analysis of miRNA and 
mRNA data 

To identify mRNAs that were potentially regulated by DE miRNAs in 
MIMs, we first sought to calculate Pearson correlations using the 
Apache commons Java statistics library between all DE miRNA 
expression in reads per million (rpm) and all mRNA expression 
(TMM normalized read counts) over the time course. A correlation 
matrix was constructed consisting of 26 (DE miRNAs) X 24,616 (all 
mRNAs) correlation coefficients. The resulting correlation matrix was 
then filtered to remove nonsignificant correlations (critical value for 
Pearson correlation for this matrix is r = —0.3116) and those inverse 
correlations that were not supported by miRanda-predicted miRNA- 
target pairs (miRanda v3.3a). The mRNAs predicted to be targeted 
{i.e., had a significant anticorrelation relationship in the expression of 
the mRNA and the miRNA, plus a predicted seed target) by either 
upregulated or downregulated miRNAs were then selected for path- 
way analysis. Two-dimensional cluster analysis and visualization, us- 
ing R version 2.15.3 hclust and heatmap.plus packages, were then 
applied to the filtered correlation matrix. 

Pathway analysis of predicted miRNA target genes 

Target genes of DE miRNAs were submitted to InnateDB (Lynn et al 
2008) for pathway analysis. Genes were submitted in two groups, 
those that were targets of upregulated miRNAs and those that were 
targets of downregulated miRNAs. Significant pathways were calcu- 
lated based on hypergeometric analysis; pathways of interest were 
defined as having a Benjamini and Hochberg corrected P value 
of < 0.05 (Benjamini and Hochberg 1995). 

Novel miRNA discovery 

Using the software package miRDeep2 vO.0.5 (Mackowiak 2011), we 
examined whether milk/blood monocytes encoded for miRNAs not 
yet annotated in the bovine genome. We further parsed these data 
using a number of different parameters to identify those novel 
miRNAs that have the highest likelihood of being true positives as 
described previously (Lawless et al 2013). Specifically, we identified 
those predictions when the mature and star stands were expressed 
with a minimum of five reads each; when miRDeep2 predicted that 
the miRNA had >90% probability of being a true positive; when the 
hairpin structure had a significant Randfold P value; and when the 
novel miRNA was independently predicted in two or more different 
miRNAseq samples. 
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RESULTS 

Kinetics of 5. uberis infection in vivo 

To investigate the host monocyte transcriptional and posttranscrip- 
tional responses to a mastitis-causing pathogen in vivo, five Holstein 
Friesian animals were infected via the teat canal with approximately 
500 CFU of Streptococcus uberis 0140 in 10 ml saline. Five control 
animals were inoculated with saline only. Blood and milk samples were 
taken at 0, 12, 24, 36, and 48 hr postinfection (hpi) and CD 14+ 
monocytes were isolated by FACS (Figure 1). The infection was mon- 
itored using recorded milk bacterial counts (CFU/ml) and somatic cell 
counts (per ml) at each of the five time points for each animal (control 
and infected). On average, bacterial counts peaked in the infected 
animals at 24 hpi, whereas no change was observed in the uninfected 
controls. There was, however, significant heterogeneity among the 
CFU data for each infected animal in terms of the magnitude and 
the timing of the response (Figure 2A and Table S3). One infected 
animal (TI3, which had the highest CFU/ml data) peaked at 24 hpi, 
two others (Til and TI4) peaked at 36 hpi, and CFU data for one 
animal were still climbing at 48 hpi (TI5). Additionally, one infected 
animal was observed to have only a very modest increase in bacterial 
counts (TI2). Somatic cell count (SCC) data also confirmed the pres- 
ence of the infection in the challenged animals and not in the controls. 
SCC was observed to increase at 24 hpi and by 36 hpi was, on average, 
> 900,000 cells/ml in infected animals. In comparison, the average 
SCC in control animals at 36 hpi was < 52,000 cells/ml (Figure 2B 
and Table S3). An SCC reading >200,000 cells/ml is generally con- 
sidered diagnostic of mastitis (Dufour and Dohoo 2013). Again, there 
was heterogeneity in the SCC response in the infected animals. In- 
terestingly, the infected animal that was observed to have only a modest 
increase in CFU/ml (TI2) had a relatively robust SCC response. 

Profiling mRNA expression in blood-isolated and milk- 
isolated 0014"^ monocytes 

NGS approach was applied to monitor the transcriptional (mRNAseq) 
and posttranscriptional (miRNAseq) changes that occurred in blood- 
isolated and milk- isolated CD 14+ monocytes in the infected and con- 
trol animals (Figure 1). Sequencing of 100 mRNA lUumina libraries 
{i.e., 50 blood and 50 milk monocyte mRNA libraries) yielded more 
than four billion sequence reads. More than three billion reads of these 
mapped uniquely to the Bos taurus UMD 3.1 genome (Table S2). 



The average correlation coefficient of mRNA normalized read counts 
between samples at each time point was 0.95 for control samples and 
0.92 for infected samples, indicating very high reproducibility of the 
data among replicates (Table S4). 

Hierarchical clustering of normalized mRNA read counts from 
MIMs revealed that the control and infected animals clearly separated 
at 36 and 48 hpi, except for the one infected sample (TI2) that had very 
low bacterial counts and likely did not develop a full infection (Figure 
2, C-F). This sample was subsequently excluded from differential gene 
expression analysis. Hierarchical clustering of the normalized mRNA 
read counts of genes that were DE in blood-isolated monocytes (BIMs) 
revealed that only three of the infected animals (Til, TI3, and TI4) 
separated from uninfected controls at 36 and 48 hpi (Figure SI). These 
animals also had the highest SCC data at these time points. 

We utilized the EdgeR statistical package (Robinson et al. 2010) to 
determine which mRNAs were significantly DE in MIMs and BIMs in 
response to S. uberis. In MIMs, there were 4, 36, 1774, and 1532 upre- 
gulated genes at 12, 24, 36, and 48 hpi, respectively. The majority (1254) 
of genes upregulated at 48 hpi were also upregulated at 36 hpi. Addi- 
tionally, there were 5, 2, 1518, and 995 downregulated genes in MIMs at 
those time points. Of the 995 downregulated genes at 48 hpi, 80% were 
also downregulated at 36 hpi. Overall, 2056 different genes were upre- 
gulated and 1721 different genes were downregulated for at least one 
time point in MIMs in response to S. uberis infection (Table S5). 

Traditionally, mastitis has been thought of as a local bacterial 
infection with a robust inflammatory response. In BIMs, however, we 
observed a much more subtle but still quite significant response to S. 
uberis infection. Only 10 genes were upregulated in BIMs at 36 hpi, but 
this increased to 83 genes by 48 hpi (Table S5). Nine of the 10 36 hpi 
genes were also upregulated at 48 hpi. Additionally, 3, 4, 26, and 39 
genes were downregulated in BIMs at 12, 24, 36, and 48 hpi, respectively. 

Pathway analysis reveals the suppression of metabolic 
pathways and the upregulation of inflammatory 
pathways in response to 5. uberis infection 

Pathway analysis of upregulated and downregulated genes at each 
time point was undertaken using GOseq (Young et al. 2010) with 
pathway annotation imported from the KEGG database (Kanehisa 
et al. 2007) to identify which pathways were statistically overrepre- 
sented among DE genes in MIMs and BIMs. Two manually curated 
pathways (interferon signaling pathway and the inflammasome 
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Figure 1 Overview of the experimental 
design. Five Holstein Friesian animals 
were infected via the teat canal of the 
right front quarter with approximately 
500 CFU of a mastitis-causing pathogen, 
Streptococcus uberis 0140, in 10 ml sa- 
line. Five control animals were inocu- 
lated with saline only. Milk and blood 
samples were obtained from each animal 
at 0, 12, 24, 36, and 48 hr after infection 
(or mock infection). Milk-derived and 
blood-derived CD14+ monocytes were 
isolated by fluorescence-activated cell 
sorting (FACS) from each sample. mRNA 
and miRNA were extracted and 200 
lllumina-compatible libraries were pre- 
pared for sequencing on a Hiseq 2000 
machine. More than four billion mRNA 
reads and more than two billion miRNA 
reads were sequenced in total. 
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Figure 2 The response to infection. (A) The infection was monitored using recorded milk bacterial counts (CFU/ml) and (B) somatic cell counts 
(per ml) at each of the five time points (0, 1 2, 24, 36, and 48 hpi) for each animal (control and infected). Significant heterogeneity was observed in 
the CFU data among the infected animals. One infected animal (TI2) was observed to have only a very modest increase in bacterial counts. (C-F) 
Hierarchical clustering of the top 500 most variable probes (more red the colour the more highly expressed the gene is) in milk-isolated 
monocytes (MIMs) at 12, 24, 36, and 48 hpi, respectively, revealed that infected animals separated from control animals in their gene expression 
response at 36 and 48 hpi except for animal TI2. (G-J) InnateDB network analysis of genes that were differentially expressed at 1 2, 24, 36, and 48 
hpi, respectively (red = upregulated: green = downregulated). The networks were visualized in Cytoscape analysis of genes that were 
differentially expressed at 12, 24, 36, and 48 hpi, respectively. The networks were visualized in Cytoscape. 



pathway) were also included (see Materials and Methods). No sig- 
nificant pathways were identified among either the BIM or the MIM 
DE genes at 12 or 24 hpi. At 36 and 48 hpi, however, more than 20 
different pathways were identified as being statistically overrepre- 
sented (Figure 3 and Table S6). In MIMs, downregulated genes were 
predominantly associated with metabolic pathways (Figure S2), such 
as fatty acid and amino acid metabolism, the citric acid (TCA) cycle 
and glutathione metaboUsm, as well as DNA replication and repair 
pathways and the cell cycle. Downregulated pathways were largely 
similar between 36 hpi and 48 hpi, although fewer pathways were 
significant at 48 hpi. Upregulated genes, however, were primarily 
associated with well-known pattern recognition receptor pathways 



(Figure 4), including the Toll-like receptor pathway [e.g., TLR2, 
TLR4, CD14, MYD88, TIRAP, and IL-1 receptor-associated kinase 1 
(IRAKI) all upregulated], the NOD-like receptor pathway [e.g., NODI, 
NOD2, NLRP3 (NALP3), NLRC4 (IPAF), NAIP (NAIP5) upregulated], 
and the RIG-I-like receptor pathway [e.g., DDX58 (RIG-I), IFIHl 
(MDA5), CYLD, DHX58 (LGP2), DDX3X, TRIM25], and interferon 
signaling and cytokine and chemokine signaling pathways. Upregulated 
inflammatory cytokine and chemokine genes included the genes encod- 
ing TNF, IL-IA, IL-IB, IL-6, IL-8, IL-12A and IL-12B, IL-17B and 
IL-17C, IL-18, IL-23A, IL-27, CCL3 (MlPla), CCL4 (MIPip), 
CCL5 (RANTES), CCL8 (MCP-2), and CCL20 (MIP3A). The genes 
encoding TNF, IL-IB, IL-6, IL-12, and CCL20 were more flian 10-fold 
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Figure 3 Pathways that were statistically 
overrepresented among genes. (A) Upre- 
gulated at 36 hpi, (B) upregulated at 48 
hpi, (C) downregulated at 36 hpi, and (D) 
downregulated at 48 hpi in MIMs iso- 
lated from S. uberis- infected animals. 
Upregulated genes were significantly en- 
riched for roles in inflammatory and other 
innate immune pathways, whereas downre- 
gulated genes were significantly associated 
with metabolic pathways. 



upregulated at 36 hpi. All upregulated pathways that were significant at 
36 hpi were still significant at 48 hpi, with five additional upregulated 
pathways being significant only at 48 hpi. These pathways were primar- 
ily related to leukocyte migration and phagocytosis. 

In BIMs, only two pathways were statistically overrepresented among 
48 hpi upregulated genes, interferon signaling and cytokine-cytokine 
receptor interaction (Table S6). No pathways were significant at the other 
time points. Among downregulated genes, there were also few overrep- 
resented pathways. Those pathways that were significant were primarily 
related to the complement and focal adhesion pathways. 

Network analysis of differentially expressed genes 

InnateDB was used to generate a network of experimentally supported 
molecular interactions that have been annotated to occur directly 
between the DE genes and their encoded products. Gene expression 
data from each of the four postinfection time points were then overlaid 
on this network and the network was visualized using Cytoscape 2.8.2 
(Shannon et al. 2003) (Figure 2, G J). The network consisted of 2185 
nodes (representing DE genes and their encoded products) and 10,786 
edges (representing annotated molecular interactions) between them. 



The network was then analyzed using cytoHubba (Lin et al. 2008) 
to identify network hubs and bottlenecks that may represent the key 
regulatory nodes in the networks. Using the "Degree" algorithm, the 
top 20 hubs {i.e., genes/proteins that are highly connected to other DE 
genes) in each network were identified (Table S7). InnateDB Gene 
Ontology analysis revealed that the top 20 hubs were highly enriched 
for roles in innate immunity (FDR <1.7e~^) and transcriptional reg- 
ulation (FDR <1.5e~^). Many of the top 20 hubs were well-known 
transcriptional regulators of innate immunity, including JUN, NFKBl, 
RELA, STATl, STAT3, EP300, and CREBBP. These transcriptional 
regulators were located in the most densely connected portion of the 
network and share many connections (Figure 5A). One interpretation 
of this is that transcriptional regulation by several of these transcrip- 
tion factors is required for differential gene expression in response to 
the infection. 

Bottlenecks are network nodes that are the key connector proteins 
in a network and have many "shortest paths" going through them, 
similar to bridges or tunnels on a highway map (Yu et al. 2007). 
Seventeen of the top 20 hubs were also identified as network bottle- 
necks. Beta-actin, the transcriptional regulator, FOS (API), and ISG15 
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Figure 4 Differential gene expression in 
key innate immune pathways in MIMs 
from S. ubens- infected animals. (A-C) 
Differentially expressed (DE) genes high- 
lighted on KEGG Toll-like receptor (TLR), 
NOD-like receptor (NLR), and RIG-I signal- 
ing pathway diagrams (red = upregulated; 
green = downregulated). (D-F) RPKM 
values for each of the DE genes in the 
TLR, NLR, and RIG-I pathways, respectively 
(red = infected; blue = control). (C) Fold 
change values for each of the DE genes 
in the TLR, NLR, and RIG-I pathways, 
respectively (red = upregulated; green = 
downregulated). 



TLR Pathway - 36hpi Fold change 



NLR Pathway - 36hpi Fold change 



RIG-I Pathway - 36hpi Fold change 
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Figure 5 Network analysis of differentially expressed (DE) genes. InnateDB was used to generate a network of experimentally supported 
molecular interactions that have been annotated to occur directly between the DE genes and their encoded products. (A) The top 20 CytoHubba- 
predicted network hubs and their interactors. Overlain on the network is gene expression data from MIMs at 36 hpi (red = upregulated; green = 
downregulated). The top 20 hubs were highly enriched for roles in innate immunity. (B) The high-scoring DE subnetwork identified in the larger 
network using jActive Modules. This module was highly enriched for several innate immune relevant pathways. 



were additionally identified in the top 20 bottlenecks. ISG15 has been 
shown to act as a negative regulator of both the NF-kB and RIG-I 
signaling pathways (Arnaud et al 2011; Minakawa et al 2008). 

Identifying hubs solely based on their degree can identify nodes 
that, in general, have been annotated to have many interactions. This 
fails to address whether the number of connections in a particular 
network of interest is more or less than is expected (given the number 
of known interactions for that node in the database and the size of the 
network). To address this issue, we have developed the Contextual 
Hub Analysis Tool (CHAT) (H., L. Wiencko, K., Bryan, A., Noyelle- 
Depierre & D., J. Lynn., unpublished results ), a network analysis tool 
that identifies nodes in a network that are more highly connected to 
contextually relevant nodes (in this case DE nodes) than is expected 
by chance (Table S8). Applying this method to our network identified 
that several of the hub nodes (e.^., KIAAOlOl) did not have more 
connections to DE genes than expected by chance. This means that 
these genes, although highly connected, are less likely to be function- 
ally important in our network. However, all of the transcriptional 
regulators of innate immunity that were identified as hubs in the 
analysis above (JUN, NFKBl, RELA, STATl, STAT3, EP300, 
CREBBP), were also identified by CHAT to be significantly more 
connected to 36 hpi DE genes than expected by chance. CHAT also 
identified several other known innate immunity transcriptional regu- 
lators in the top 20 "contextual' hubs, including REL, IRFl, and IRF9. 
Fifteen of the top 20 contextual hubs are annotated by InnateDB as 
having a role in innate immunity (FDR <1.16E~^^). Although the 
ranking of the top contextual hubs changed from 36 to 48 hpi, the 
most significant hubs remained the same. 

The network was also analyzed using the jActiveModules plugin 
(Ideker et al 2002) in Cytoscape 2.8.2 (Shannon et al 2003) to identify 
high-scoring DE subnetworks. This type of analysis can aid in the 
identification of functionally relevant groups of DE genes that may 
be acting in concert. A single highly connected component (>10 
nodes) was identified consisting of 278 nodes and 1585 interactions. 



This module consisted of several of the transcriptional hubs (CREBBP, 
EGRl, EP300, NFKBl, RELA) identified in the analysis above and 
their interactors (Figure 5B). Pathway analysis of genes in the module 
revealed that the module was statistically enriched for many of the 
same pathways that were identified in the analysis of all upregulated 
genes including Jak-STAT signaling, the TLR, NLR, and RIG-I path- 
ways, apoptosis, and chemokine and cytokine signaling (Table S9). 

Profiling miRNA expression in blood-isolated and milk- 
isolated 0014"^ monocytes 

Sequencing of 100 miRNA lUumina libraries yielded more than one 
biUion reads for both the MIM and BIM samples. Following a pipeline 
of quality filtering and adaptor trimming, a total of 312 and 492 
million reads from MIMs and BIMs, respectively, mapped uniquely to 
the Bos taurus UMD 3.1 genome (Table S2). Uniquely aligning reads 
were then assigned to known mRNAs/miRNAs using HTseq based on 
miRBase vl9 annotation of the bovine genome to generate read 
counts per mature miRNA in each sample (Flicek et al 2012; Hubbard 
et al 2009). On average, 79% of MIM reads and 80% of BIM reads 
uniquely mapped to known miRNAs (Figure S3). 

RNAseq profiling revealed that the miRNAome of MIMs and BIMs 
were broadly similar and exhibited a range of miRNA expression that 
has been observed in many other cell types; 297 and 282 miRNAs were 
expressed at a threshold of >1 rpm in MIMs and BIMs, respectively 
(Table SIO). Of these, 136 and 116, respectively, were expressed at 
a level >100 rpm, a level of expression that has been shown to be 
associated with functional miRNAs (Mullokandov et al 2012). 

Multiple miRNAs are differentially expressed in 
response to S. uberis infection 

The EdgeR statistical package was utilized to determine which 
miRNAs were significantly DE in response to S. uberis infection at 
12, 24, 36, and 48 hpi. To address any normalization issues, we 
retained only those miRNAs that were robust to the normalization 
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procedure used (Lawless et al. 2013). Additionally, only miRNAs that 
had an average expression of >10 rpm across samples were included 
for further analysis. In MIMs, we identified that 26 unique miRNAs 
were DE. Twelve of these miRNAs were DE across more than one 
time point (Table 1). Hierarchical clustering of the normalized read 
counts for MIM DE miRNAs revealed that the control and infected 
animals clearly separated at 36 hpi (Figure 6). Very few miRNAs were 
DE in blood monocytes. We found three in total, one at 24 hpi and 
three at 48 hpi (one was DE at both time points) (Table 1). 

Many of the miRNAs identified as being DE have been shown to 
have a role in immunity in other species. miR-223, for example, which 
was upregulated at 36 hpi in MIMs, has a multifactorial role in 
neutrophils, regulating their proliferation, activation, and granulopoiesis 
(Chen et al 2004; Fazi et al 2005). Interestingly, in RAW264.7 cells 
challenged with lipopolysaccharide (LPS), miR-223 has been reported to 
be downregulated, allowing the upregulation of signal transducer and 



activator of transcription 3 (STAT3), which promotes proinflammatory 
IL-6 and IL-ip transcription. We have previously commented on the 
differences in the miRNA response to LPS and S. uheris, a Gram- 
positive bacterium (Lawless et al 2013). 

Other DE miRNAs in our study that have a demonstrated role in 
immunity and infection in other species include let-7e, which was 
downregulated in MIMs after S. uheris infection. Let-7e has been 
shown to regulate caspase 3 and caspase 7 in human monocyte- 
derived macrophages infected with Mycobacterium avium hominissuis 
(Sharbati et al 2011). Let-7e, as well as several other DE miRNAs in 
our study (bta-miR-200c, bta-miR-210, and bta-miR-193a), have also 
previously been identified as DE in bovine mammary epithelial cells 
stimulated with S. uheris in vitro (Lawless et al 2013). 

Another DE miRNA with a role in immune regulation is miR-150, 
which we observed to be downregulated at both 36 and 48 hpi in MIMs. 
miR-150 targets MyD88 (upregulated at both 36 and 48 hpi in MIMs), 



W Table 1 Fold changes and false discovery rates of DE miRNAs at 12, 24, 36, and 48 hr postinfection in milk-isolated and blood-isolated 
monocytes 
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Figure 6 A heatmap of the normalized read counts of miRNAs that were 
identified as differentially expressed in MIMs at 36 hpi. Hierarchical 
clustering revealed the separation of the infected and control animals 
based on this miRNA expression data (more red the colour the more 
highly expressed the miRNA is). Note that one infected animal (TI2) is 
not included in this analysis because it did not appear to respond to the 
infection as measured by CFU or mRNA expression data (see Figure 2). 



a key regulator of TLR signaling (Ghorpade et al 2013). miR-lSO has also 
been shown to target CXCR4 (Rolland-Tumer et al 2013; Tano et al 
2011), which was upregulated two-fold at 36 and 48 hpi in MIMs. Finally, 
one of the miRNAs that was identified as downregulated in BIMs, miR- 
146b, has been shown to target TNF receptor-associated factor 6 and 
IRAKI genes in THP-1 cells stimulated with LPS (Taganov et al 2006). 

Predicted targets of downregulated but not 
upregulated miRNAs are highly enriched for roles in 
innate immunity 

To identify the potential mRNA targets of DE miRNAs in MIMs 
isolated from infected animals, we identified those mRNAs whose 
expression was significantly negatively correlated with miRNA 
expression. These predictions were further refined by removing those 
miRNA target-predicted relationships that were not supported by 
a predicted seed region in the 3' UTR of the correlated mRNA (Figure 
7 and Table SU). Analysis of the predicted target genes using 
InnateDB (www.innatedb.com) (Lynn et al 2008) revealed that the 
predicted targets of downregulated but not upregulated miRNAs were 
highly enriched for roles in innate immunity (FDR <3.2E~^). More 
specifically, pathway analysis revealed that downregulated miRNAs 
were predicted to preferentially target key pathogen recognition re- 
ceptor signaling pathways, including the TLR, NLR, and RIG-I sig- 
naling pathways (Table SI 2). Given that these pathways were also 
identified in the mRNAseq data as being among the top upregulated 
pathways, this finding strongly suggests that miRNAs are key regu- 
lators of innate immune pathways that drive the host inflammatory 
response during mastitis. 



In contrast, the predicted targets of upregulated miRNAs were 
enriched for roles in metabolism (FDR = 0.01). Pathway analysis of 
the mRNAseq data had already highlighted the downregulation of 
metabolic pathways in response to S. uheris infection (see above). 
These results suggest that miRNAs may also be key regulators of 
the transcriptional suppression of metabolic pathways during mastitis. 

Novel miRNA discovery 

The miRNA sequencing data were also mined to determine if milk or 
blood monocytes expressed potentially novel miRNAs. Previously, we 
identified 19 novel miRNAs in bovine mammary epithelial cells using 
miRDeep2 (Lawless et al 2013). Applying the same approach as our 
previous study, we identified a further 20 high- confidence putatively 
novel bovine miRNAs that were independently predicted in multiple 
MIMs/BIMs miRNAseq data (Table S13). 

Searching the miRBase database (v20) using BLAST identified that 
eight of the novel miRNAs had close homology to other miR-2284 
family members. miRBase currently lists this miRNA family as having 
102 members, yet virtually no data exist regarding what function such 
a large group of constitutively expressed miRNAs may have. Other 
novel miRNAs discovered in this study included a homolog of hsa- 
miR-3680-3p, a miRNA identified in human periphery blood (Vaz 
et al 2010). The discovery of these miRNA further adds to the data- 
base of bovine miRNAs. 

DISCUSSION 

Infectious disease is a serious threat not only directly to human health 
but also to animal health, for which it is associated with substantial 
annual economic losses, public confidence issues, and food security 
concerns. Currently, there is a significant gap in the understanding of 
the molecular and genetic mechanisms that underpin susceptibility to 
infectious disease in humans and in animals. An important reason for 
this is the fact that disease susceptibility is a multifactorial complex 
phenotype, which is not the result of single genes acting in isolation 
but rather is attributable to perturbation at a network or systems level 
(Barabasi et al 2011). Such networks are regulated at multiple differ- 
ent levels {e.g., genetic, transcriptional, posttranscriptional) and, as 
such, a multiomic integrative biology approach is needed to under- 
stand them. 

Here, we report NGS approach coupled with advanced network 
and pathway biology methods to simultaneously profile the mRNA 
and miRNA networks that are differentially regulated in vivo in blood- 
isolated and milk- isolated CD 14+ monocytes during infection with 
a bovine mastitis pathogen, S. uheris. Bovine mastitis is an inflamma- 
tion-driven disease of the bovine mammary gland, which costs the 
global dairy industry billions of dollars per year (Jones and Bailey 
2009; Wells et al 1998). Profiling genome- wide changes in mRNA 
expression in MIMs and BIMs using RNAseq, we observed more than 
3500 genes to be statistically altered in their expression in response to 
S. uheris challenge. Notably, this RNAseq approach identified approx- 
imately 1000 more DE genes than had been previously reported in 
a microarray-based analysis of RNA expression in mammary tissue of 
cattie infected with the same pathogen (Jensen et al 2013; Moyes et al 
2009; Swanson et al 2009). As expected, given that mastitis is a rela- 
tively localized inflammatory disease in the mammary gland, the ma- 
jority of DE genes were identified in MIMs, which are recruited to the 
site of infection. This influx of immune cells (including monocytes) to 
the site of infection was observed in recorded SCC data only in the 
infected animals. However, we also observed a small but significant 
transcriptional response in BIMs to S. uheris infection that was pri- 
marily associated with an interferon and cytokine signaling signature. 
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Figure 7 Visualization of miRNA targets. (A) Two-dimensional cluster analysis and heatmap visualization of the miRNA/mRNA expression 
correlation data matrix after filtering by predicted miRNA-targets and correlation significance. The primary split in the upper hierarchical 
dendrogram largely aligns with the upregulated miRNA (red) and downregulated miRNA (green). miRNA/mRNA target expression correlations are 
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Previous studies have also shown more systemic changes in gene 
expression in neighboring uninfected mammary glands, the liver, 
and in the blood (Blum et al 2000; Jensen et al 2013; Jiang et al 
2008; Mitterhuemer et al 2010). 

The predominant signature associated with upregulated mRNAs in 
MIMs from S. w^er/s-infected animals was the strong transcriptional 
activation of innate immune and inflammatory gene expression. In 
particular, we noted the transcriptional activation of key pattern rec- 
ognition pathways, including the TLR, NLR, and RIG-I pathways, 
which likely drive the observed proinflammatory response. The in- 
volvement of TLR signaling (particularly TLR2 and TLR4) in the host 
response to mastitis is well-documented (Buitenhuis et al. 2011; Ma 
et al 2011; Mitterhuemer et al 2010; Porcherie et al 2012; Whelehan 
et a/. 201 1); however, less is known about the involvement of the NLR 
and RIG-I pathways (Moyes et al 2009). Interestingly, the RIG-I 
pathway is classically associated with viral RNA recognition; however, 
recent findings suggest that RIG-I can also recognize nucleic acids 
released by invasive bacteria and trigger IFN-P and inflammasome 
activation (Abdullah et al 2012). These findings concur well with the 
observed interferon and inflammasome activation transcriptional sig- 
natures observed in our study. 

Several previous in vitro studies have strongly suggested roles for 
miRNAs in regulatiag bovine immunity (Dilda et al 2011; Lawless et al 
2013); however, none of these have globally profiled the miRNA re- 
sponse to infection in vivo. In this study, we have also used NGS 
sequencing to profile miRNA expression in MIMs and BIMs after S. 
uheris infection. Twenty-six miRNAs were identified as DE in MIMs 
and three were identified in BIMs. Several of these have been previously 
described as targeting immune or inflammatory regulators in other 
species. Of particular interest is our finding that downregulated but 
not upregulated miRNAs in MIMs are predicted to preferentially target 
genes involved in innate immunity and inflammation. Furthermore, the 
TLR, NLR, and RIG-I pathways discussed above were aU preferentially 
predicted to be targeted. This strongly suggests that the transcriptional 
suppression of these miRNAs enables the activation and amplification 
of the proinflammatory response. Further supporting this conclusion is 
the fact that several of the DE miRNAs in our study have been validated 
to target genes in these pathways in other species. miR-149, for example, 
which was downregulated in MIMs after S. uheris infection, has been 
shown to target mouse CD 14 and IRAKI (Chi et al 2009), key signaling 
proteins in the TLR pathway. Both CD 14 and IRAKI are also predicted 
to be targets of bta-miR-149 in our study. 

The other predominant transcriptional signature that we found in 
MIMs after S. uheris infection was the widespread repression of a num- 
ber of metabolic processes (>150 KEGG-annotated metabolism genes 
are downregulated at 36 hpi). Interestingly, we found that upregulated 
miRNAs were predicted to preferentially target genes involved in 
metabolism, suggesting that miRNAs, which are upregulated in re- 
sponse to S. uheris infection, may contribute to the transcriptional 
suppression of metabolic pathways. This signature of metabolic gene 
transcriptional suppression may appear initially paradoxical in light of 
the fact that production of inflammatory cytokines is expected to 



require substantial energy consumption. It is now becoming widely 
appreciated that activated macrophages undergo the Warburg effect, 
switching their metabolism from oxidative phosphorylation to glycol- 
ysis (McGettrick and O'Neill 2013). This metabolic switch has recenfly 
been investigated using a combined metabolomics and microarray 
approach (Tannahill et al 2013) and has revealed the upregulation 
of a number of genes involved in glycolysis in bone marrow-derived 
macrophages challenged with LPS [e.g., solute carrier family 2 (facilitated 
glucose transporter), member 1 (SLC2A1/GLUT1), hexokinase 
3 (HK3), fructose-2,6-biphosphatase 3 (PFKFB3)] and the downregula- 
tion of several key genes encoding enzymes in the TCA cycle [e.g., 
malate dehydrogenase 1 (MDHl) and isocitrate dehydrogenase 2 
(IDH2)]. Our transcriptional data are consistent with the data presented 
in fliis article {e.g., SLC2A1; HK2 and HK3 and PFKFB4 are all 
upregulated and MDHl and IDHl are downregulated in MIMs at 36 
hpi) and suggest that although there is a broad signature of transcrip- 
tional suppression of metabolism, these cells are likely to be highly 
glycolytically active. Several other genes encoding enzymes in the TCA 
cycle (which was statistically overrepresented among downregulated 
genes) were also transcriptionally repressed in MIMs at 36 hpi, in- 
cluding the following: dihydrolipoamide dehydrogenase; fiamarate 
hydratase; IDH3B; pyruvate dehydrogenase beta; MDH2; succinate 
dehydrogenase complex, subunit B; and succinate-CoA ligase, alpha 
subunit. As has also recenfly been shown in bone marrow-derived 
macrophages, LPS strongly increases levels of succinate, a TCA cycle 
intermediate (Tannahill et al 2013). Succinate acts as an inflammatory 
signal in macrophages inducing ILIB through the transcription factor 
HIF-la, both of which are transcriptionally activated in MIMs 36 hpi 
{ILlB is 10-fold upregulated). Another metabolic pathway that is signif- 
icanfly transcriptionally repressed in MIMs at 36 and 48 hpi is the KEGG 
Valine, leucine, and isoleucine degradation pathway. Valine, leucine, and 
isoleucine are branch-chain amino acids that are converted into 
Acyl-CoA derivatives. These are converted into either acetyl-CoA or 
succinyl-CoA and enter the TCA cycle (Sears et al 2009). Most of the 
other significantiy downregulated pathways, including fatty acid metab- 
olism, propanoate metabolism, butanoate metabolism, tryptophan 
metabolism, beta- alanine metabolism, lysine degradation, and 
glyoxylate and dicarboxylate metabolism, also result in the pro- 
duction of acetyl- Co A and succinyl-CoA that enter the TCA cycle. 
A similar pattern of expression leading to the transcriptional 
downregulation of the TCA cycle and alternative pathways in- 
volved in producing TCA cycle components has also been reported 
in other infection models (Chin et al 2010). The transcriptional 
repression of these pathways is therefore also consistent with 
a switch in metabolism from oxidative phosphorylation to glycol- 
ysis during the proinflammatory response. 

Another pathway that is significantiy downregulated is the KEGG 
primary bile acid biosynthesis pathway. Despite the potentially 
misleading name, this pathway primarily consists of the reactions 
involved in cholesterol metabolism. The downregulation of genes 
involved in cholesterol metabolism has also been reported in 
monocytes isolated from HIV^ individuals (Feeney et al 2013) and 



colored based on increasing significance (light blue to blue) with nontarget/nonsignificant correlations masked (gray). (B) A network 
representation of the predicted targets (circular nodes) of downregulated miRNAs (green arrows). Red circles = upregulated in MIM mRNA 
expression data at 36 and/or 48 hpi. Larger circular nodes represent those genes that have been annotated by InnateDB to have a role in innate 
immunity. (C) A network representation of the predicted targets (circular nodes) of upregulated miRNAs (red arrows). Green circles = 
downregulated in MIM mRNA expression data at 36 and/or 48 hpi. Note that the predicted targets of downregulated but not upregulated 
miRNAs are highly enriched for roles in innate immunity. 
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a number of bacterial infections (Dushkin 2012). The accumulation of 
cholesterol in monocytes and macrophages leads to the formation of 
foam cells, which in humans are associated with the inflammatory 
disease, atherosclerosis (Ross 1999). The downregulation of genes in- 
volved in cholesterol and fatty acid metabolism is likely driven in part 
by the transcriptional suppression of the PPAR-7 transcription factor 
in MIMs at 36 hpi and the downregulation of PPAR-a at 48 hpi, both 
of which are key transcriptional regulators of these pathways 
(Dushkin 2012). Interestingly, the PPARs also have a role in the 
regulation of inflammation, in which their suppression is required 
to induce inflammatory gene expression (Bensinger and Tontonoz 
2008). Of further note is that one of the liver X receptors, LXR-P 
(NR1H2), which together with the PPARs is a key regulator of in- 
flammation and lipid metabolism (Bensinger and Tontonoz 2008), is 
upregulated in MIMs at 36 and 48 hpi. LXRs are also known to 
antagonize inflammatory gene expression, so it is somewhat surprising 
to find LXR-P to be upregulated. This may reflect the fact that balance 
is needed to avoid excessive inflammation or that LXRs and PPARs do 
not completely overlap in the genes that they regulate. 

An additional link between metabolism and inflammation that is 
currently undergoing intensive investigation is the role of NAD+, 
sirtuins (SIRTs), and AMP-dependent protein kinase (AMPK) in 
suppressing inflammation (McGettrick and O'Neill 2013). The acti- 
vation of TLR4, which along with TLR2 and TLR9 is transcriptionally 
upregulated in MIMs at 36 hpi, has been shown to induce NAM 
phosphoribosyltransferase (NAMPT; upregulated in MIMs at 
36hpi), which in turn activates SIRTl (upregulated in MIMs at 36 
hpi) via NAD+. SIRTl limits inflammation by repressing RELA 
transcription factor activity, a key transcriptional hub identified in 
MIMs. SIRTl also activates AMPK (upregulated in MIMs at 36 hpi), 
a central regulator of energy metaboUsm. Activation of AMPK has 
been shown to decrease NF-kB activity and TNF-a production in 
macrophages stimulated with LPS, IL-12 production in DCs, and 
HIF-la. These data suggest that at 36 hpi, the brakes are starting 
to be applied to limit the inflammatory response to S. uheris in- 
fection via SIRTl and AMPK. The affect of this break is apparent 
at 48 hpi, where the genes encoding TNF-a, IL-IB, IL-12, and HIF- 
la are all downregulated in comparison to 36 hpi. The limiting of 
inflammation at this stage makes sense based on the bacterial count 
data, which show that bacterial CFU/ml are declining, suggesting 
that the infection is being resolved. Interestingly, AMPK activity also 
inhibits both the cholesterol and fatty acid metabolic pathways that, 
as discussed above, are downregulated in MIMs after S. uheris in- 
fection. This suppression of fatty acid metaboUsm has been shown to 
be beneficial to the host during a viral infection (Moser et al. 2012). 

Finally, we also observed the statistically significant over- 
representation of downregulated genes annotated in the KEGG 
glutathione metabolism pathway. Glutathione has an important role 
in innate and adaptive immunity and has been shown to confer 
protection against microbial, viral, and parasitic infections (Morris 
et al. 2013). Glutathione metabolism also plays an important role in 
macrophages in the detoxification of reactive oxygen species. The 
transcriptional suppression of this pathway in MIMs after S. uheris 
infection may be a consequence of the metabolic switch to glycolysis 
and may be detrimental to the host, leading to an excess of oxidants in 
the cells, which could drive the inflammation and tissue damage that 
are characteristic of mastitis. In patients with active tuberculosis, 
PBMC intracellular glutathione levels declined by 70%; this was cor- 
related with increased proinflammatory cytokines and enhanced bac- 
terial growth (Guerra et al. 2012). Supplementation with glutathione 
has been demonstrated to lead to the control of mycobacterial growth 



(Venketaraman et al. 2003) and also appears to have beneficial effects 
in reducing inflammation in HIV^ patients (Morris et al. 2012). This 
suggests that glutathione supplementation is a potential strategy to 
reduce the effects of inflammation in mastitis. 

miRNAs also likely play a key role in regulating the links between 
inflammation and metabolism observed in this study. Human SIRTl, 
for example, which as discussed above limits inflammation, has been 
shown to be a target of miR-34a (Yamakuchi et al. 2008). Our data are 
consistent with this relationship also existing in MIMs, where we have 
found bta-miR-34a to be downregulated at 36 and 48 hpi and SIRTl 
to be upregulated. Other examples of miRNAs that likely transcrip- 
tionally regulate metabolic pathways in MIMs include miR-451, which 
is upregulated in MIMs at 36 and 48 hpi and has been shown to target 
M025 in mouse heart tissue altering AMPK signaling (Chen et al. 
2012). miR-451 has also been shown to regulate the expression of 
several proinflammatory cytokines in mice in response to influenza 
infection (Rosenberger et al. 2012). 

Aside from providing new insight into the regulatory role 
miRNAs play in S. uheris infection in vivo, our study also provides 
the groundwork for a number of potential practical applications in 
veterinary medicine. miRNAs, for example, exhibit many proper- 
ties that have made them of significant interest as noninvasive 
biomarkers. miRNAs are abundantly and stably expressed in 
a range of accessible tissues, including serum, milk, urine, saliva, 
and semen, where they can be readily measured (Chen et al. 2010; 
Hata et al. 2010; Kosaka et al. 2010). Importantly, for a potential 
biomarker, miRNAs have high information content, and the ex- 
pression profile of small numbers of them have been shown to be 
diagnostic of disease (De Guire et al. 2013). The use of miRNAs as 
a clinical biomarker is most advanced in human cancer research. In 
2009, Prometheus Laboratories released an miRNA biomarker to 
accurately identify 25 different tumor types (Ajit 2012), and 
miRNA biomarkers are now available for early cancer prognosis 
from two other companies, Asuragen and Rosetta Genomics. NGS- 
based technologies, such as the approach used in this study, are 
empowering RNA expression profiling, including miRNAs, on a ge- 
nome-wide scale with unprecedented resolution, accuracy, and 
speed and at a relatively low cost (Schuster 2008). There is signif- 
icant potential to develop these approaches as diagnostics of in- 
fection in animals and also in humans. 
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